function diff = q_1_p_fun(m,param)

    if param.zeta==0
        diff = zeros(size(m));
    else
        C_0_n=param.psi-(param.s.*param.lambda.*param.psi-(1-param.psi).*param.zeta)./(param.zeta+param.s.*param.lambda);
        int=(1-param.alpha).*2./param.sigma^2.*(param.zeta+param.s.*param.lambda).*C_0_n./m./q_0_fun(m,param);
        diff = zeros(size(m)) -int;
    end

end